In [1]:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import plotly.express as px
import plotly.graph_objects as go
from scipy.stats import chisquare
from warnings import warn
from functools import lru_cache
from scipy.special import gammaln
from math import lgamma, log
from pgmpy.models import BayesianModel
from functools import partial
from IPython.display import Markdown, display
import networkx as nx
import pylab as plt
def printmd(string):
display(Markdown(string))
Scoring for Structure learning DAG¶
Commonly used scoring functions to measure the fit between model and data are Bayesian Dirichlet scores：
Scoring functions for learning Bayesian networks
- BDeu
- K2
- Bayesian Information Criterion
In [2]:
# ref: https://pgmpy.org/index.html
class BaseEstimator(object):
def __init__(self, data=None, state_names=None, complete_samples_only=True):
self.data = data
if self.data is not None:
self.complete_samples_only = complete_samples_only
self.variables = list(data.columns.values)
if not isinstance(state_names, dict):
self.state_names = {
var: self._collect_state_names(var) for var in self.variables
}
else:
self.state_names = dict()
for var in self.variables:
if var in state_names:
if not set(self._collect_state_names(var)) <= set(
state_names[var]
):
raise ValueError(
f"Data contains unexpected states for variable: {var}."
)
self.state_names[var] = state_names[var]
else:
self.state_names[var] = self._collect_state_names(var)
def _collect_state_names(self, variable):
"Return a list of states that the variable takes in the data"
states = sorted(list(self.data.loc[:, variable].dropna().unique()))
return states
def convert_args_tuple(func):
def _convert_param_to_tuples(
obj, variable, parents=tuple(), complete_samples_only=None
):
parents = tuple(parents)
return func(obj, variable, parents, complete_samples_only)
return _convert_param_to_tuples
@convert_args_tuple
@lru_cache(maxsize=2048)
def state_counts(self, variable, parents=[], complete_samples_only=None):
"""
Return counts how often each state of 'variable' occurred in the data.
If a list of parents is provided, counting is done conditionally
for each state configuration of the parents.
Returns
-------
state_counts: pandas.DataFrame
Table with state counts for 'variable'
"""
parents = list(parents)
# default for how to deal with missing data can be set in class constructor
if complete_samples_only is None:
complete_samples_only = self.complete_samples_only
# ignores either any row containing NaN, or only those where the variable or its parents is NaN
data = (
self.data.dropna()
if complete_samples_only
else self.data.dropna(subset=[variable] + parents)
)
if not parents:
# count how often each state of 'variable' occured
state_count_data = data.loc[:, variable].value_counts()
state_counts = (
state_count_data.reindex(self.state_names[variable])
.fillna(0)
.to_frame()
)
else:
parents_states = [self.state_names[parent] for parent in parents]
# count how often each state of 'variable' occured, conditional on parents' states
state_count_data = (
data.groupby([variable] + parents).size().unstack(parents)
)
if not isinstance(state_count_data.columns, pd.MultiIndex):
state_count_data.columns = pd.MultiIndex.from_arrays(
[state_count_data.columns]
)
row_index = self.state_names[variable]
column_index = pd.MultiIndex.from_product(parents_states, names=parents)
state_counts = state_count_data.reindex(
index=row_index, columns=column_index
).fillna(0)
return state_counts
In [3]:
# Importing data
GME_stock_df = pd.read_csv('/etlstage/PEE_joint/NUS_modules/CS5340/group_data/GME_stock.csv')
GME_stock_df[1:] = GME_stock_df[1:].round()
GME_stock_df = GME_stock_df[GME_stock_df.date>'2011-01-01'].round(1)
GME_stock_df['date'] = pd.to_datetime(GME_stock_df['date'])
GME_stock_df.set_index('date',inplace=True)
display(GME_stock_df.head())
# Plotting the open, close, high, low, volume and adjusted close value for the term of 12 months
fig = px.line(GME_stock_df, x=GME_stock_df.index, y=GME_stock_df.columns,
title='Plot of values for a 12 month period')
fig.update_xaxes(
dtick="M12",
tickformat="%b\n%Y")
fig.update_xaxes(rangeslider_visible=True)
fig.show()
In [4]:
# Importing data Tatanic
simple_df = pd.read_csv('/etlstage/PEE_joint/NUS_modules/CS5340/group_data/small.csv')
display(simple_df.head())
Bayesian scoring functions¶
* Compute the posteriorprobability distribution, starting from a prior probability distribution on the possible networks B, conditioned to data D, that is,P(B|D).¶
* The best network is the one that maximizes the posterior probability.¶
* For comparing the BNs, can be translated to comparsion of P(B|D) (Larger P(B|D) is better).¶
Types of Bayesian scoring functions¶
- BDeu
- K2 ## Notation ### $r_i : number of states of the finite random variable X_i$ ### $q_i : number of number of possible configurations of the parent set \prod_{X_i} of X_i $ ## =========================================================================================
1. BD(K2 scoring) -> BDe -> BDeu Score¶
BD: Bayesian Dirichlet¶
Assumptions¶
- Assumption 1. Multinomial sample: Dirichlet distribution is Multinomial distribution's conjugate prior.
- Assumption 2. Dirichlet: for all parameters from given DAG G, are follow Dirichlet distribution (repeat from 1)
- Assumption 3. Parameter independence: for given DAG G where (P(G)>0):
- global parameter independence: P($\Theta_G$|G) = $\prod_i$P($\Theta_i$|G) (features to features)
- local parameter independence: P($\Theta_i$|G) = $\prod_j$P($\Theta_{ij}$|G) $\forall i$ (sample to sample)
- Assumption 4. Parameter modularity: Given 2 DAG $G_1$ and $G_2$, if $X_i$ has same parents in $G_1$ and $G_2$ Then P($\Theta_{ij}$|$G_1$) = P($\Theta_{ij}$|$G_2$) (paramters values wont change due to different graph setup <- repeat in Assumption 3 Parameter independence)
### BD score func (log ease for computation)
- $\eta_{ijk}$ are the hyperparameters for the Dirichlet prior distributions of the parameters given the network structure
- Computationally impossible: Need to compute all hyperparameters $\eta_{ijk}$ (pseudo_counts) ### K2 score func (particular case of BD)
- from BD score: $\eta_{ij} = \sum^{r_i}_{k=1}\eta_{ijk}$
- $\Gamma(.)$ is the function Gamm where $\Gamma(c)=\int_0^\infty \mathrm{e}^{-u}\mathrm{u}^{c-1}du$
- if c is an integer: $\Gamma(c)=(c-1)!$
- let all hyperparameters $\eta_{ijk} = 1$
- log(p(G)): it is quite common to assume a uniform distribution ([1] Koller & Friedman, Probabilistic Graphical Models - Principles and Techniques, 2009 Section 18.3.6.1 (esp. page 804))
In [5]:
class K2Score(BaseEstimator):
def __init__(self, data, **kwargs):
super(K2Score, self).__init__(data, **kwargs)
def K2_score_single(self, variable, parents, is_print=False):
var_states = self.state_names[variable]
var_cardinality = len(var_states) # number of element in a (undefined set) axiomatic set theory
state_counts = self.state_counts(variable, parents)
num_parents_states = float(state_counts.shape[1])
counts = np.asarray(state_counts)
if is_print:
printmd(f'r_i = :{var_cardinality}')
printmd(f'q_i = :{num_parents_states}')
printmd('state_counts($N_{ijk}$):'+str(counts))
display(state_counts.head(5))
print('-'*100)
log_gamma_counts = np.zeros_like(counts, dtype=np.float_)
# Compute log(gamma(counts + 1))
gammaln(counts + 1, out=log_gamma_counts)
if is_print:
printmd("Let $\eta_{ijk} = 1$")
printmd("$log(\Gamma(N_{ijk} + \eta_{ijk}) / \Gamma(\eta_{ijk}))$ = $log(\Gamma(N_{ijk} + 1) / 1)$ = ")
print(f"log_gamma_counts = \n {log_gamma_counts}")
print('-'*100)
# Compute the log-gamma conditional sample size
log_gamma_conds = np.sum(counts, axis=0, dtype=np.float_)
gammaln(log_gamma_conds + var_cardinality, out=log_gamma_conds)
if is_print:
printmd("$N_{ij} = \sum_k(N_{ijk})$")
printmd("$\eta_{ij} = \sum_k(\eta_{ijk}) = k * 1$")
printmd("$log(\eta_{ij}) / \Gamma(N_{ij} + \eta_{ij}))$ = $ -log(\Gamma(N_{ij} + k))$ = ")
print(f"log_gamma_conds(-ve sign will be subtract in score computation) = \n {log_gamma_conds}")
print('-'*100)
# Compute the log P(G)
log_PG = num_parents_states * lgamma(var_cardinality)
if is_print:
print(f"num_parents_states = {num_parents_states}")
print(f"var_cardinality = {var_cardinality}")
printmd("$P(G) = no\_of\_paremts * log(\Gamma(no\_of\_states)) $")
print(f"log(P(G)) = num_parents_states * lgamma(var_cardinality) = {log_PG}")
print('-'*100)
score = (
np.sum(log_gamma_counts)
- np.sum(log_gamma_conds)
+ num_parents_states * lgamma(var_cardinality)
)
if is_print:
print(f"np.sum(log_gamma_counts) - np.sum(log_gamma_conds) + log_PG:")
return score
def score(self, model,is_print=False):
score = 0
for node in model.nodes():
if is_print:
printmd(f"## node: {node}")
printmd(f"### parent: {list(model.predecessors(node))}")
printmd(f"score: {self.K2_score_single(node, model.predecessors(node),is_print=is_print)}")
print(f"="*100)
score += self.K2_score_single(node, model.predecessors(node))
return score
In [6]:
# Importing data
simple_df = pd.read_csv('/etlstage/PEE_joint/NUS_modules/CS5340/group_data/small.csv')
display(simple_df.head())
# making a BN
BN_1 = BayesianModel([['fare','age'], ['age','passengerclass'],['passengerclass','portembarked'], ['portembarked','survived']])
BN_1_G = nx.draw(BN_1, with_labels=True,node_size=3e3,width=3)
plt.show()
K2Score(simple_df).score(BN_1,is_print=True)
Out[6]:
In [7]:
# making a BN
BN_2 = BayesianModel([['fare','passengerclass'], ['passengerclass','survived'],['age','survived'], ['portembarked','survived']])
nx.draw(BN_2, with_labels=True,node_size=3e3,width=3)
plt.show()
K2Score(simple_df).score(BN_2,is_print=True)
Out[7]:
In [8]:
display(GME_stock_df.head(5))
# making a BN
BN_GMS = BayesianModel([['open_price','adjclose_price'],
['high_price','volume'],
['low_price','volume'],
['open_price','high_price'],
['open_price','low_price'],])
nx.draw(BN_GMS, with_labels=True,node_size=3e3,width=2)
plt.show()
K2Score(GME_stock_df).score(BN_GMS,is_print=False)
Out[8]:
BDe: Bayesian Dirichlet (likelihood-equivalence ) -> BDeu¶
- from BD:
- Problem:Computationally impossible: Need to compute all hyperparameters $\eta_{ijk}$ (pseudo_counts)
Assumptions (aside of BD Assumptions 1~4)¶
- Assumption 5. Likelihood equivalence
- Definition 1. Equivalent directed acyclic graphs:
- Two directed acyclic graphs are equivalent if they can encode the same joint probability distributions.
- Assumption:
- Given 2 DAG $G_1$ and $G_2$, such that P($G_1$)>0 and P($G_2$)>0, if $G_1$ and $G_2$ are equivalence then P($\Theta_D$|$G_1$) = P($\Theta_D$|$G_2$)
- Definition 1. Equivalent directed acyclic graphs:
- Assumption 6. Structure possibility
- The skeleton of any DAG is the undirected graph resulting from ignoring the directionality of every edge.
- Complete directed acyclic graph = skeleton is complete (regardless of the direction)
- Assumption:
- For any complete directed acyclic graph G, we have that P(G) > 0
BDe computation¶
- $\eta_{ijk}$ can be solved as
- $\eta_{ijk} = N' \times p(x_{ik},w_{ij}|G_0)$
- where $p(.|G0)$ represents a probability distribution associated with a prior Bayesian network G0
- $N'$ is a parameter representing the equivalent sample size.(expresses the strength of our belief in the prior distribution)
- Computationally impossible: (pseudo_counts)
- Better than BD: (may need to discuss)
- BD: $\eta_{ijk}$
- BDe: $x_{ik},w_{ij}$
- Still Need to know $p(x_{ik},w_{ij}|G_0)$ for all i,j,k
- Better than BD: (may need to discuss)
BDeu score¶
- $p(x_{ik},w_{ij}|G_0) = \frac{1}{r_{i}q_{i}}$
- the prior network assigns a uniform probability to each configuration
- $p(x_{ik},w_{ij}|G_0)$ = $\frac{1}{r_{i}q_{i} - 0}$ $for$ $x \in [0, r_{i}q_{i}]$ else $0$
- $p(x_{ik},w_{ij}|G_0)$ = $\frac{1}{r_{i}q_{i} - 0}$ $for$ $x \in [0, r_{i}q_{i}]$ else $0$
- This score only depends on one parameter, the equivalent sample size N'!
- Since there are no generally accepted rule to determine the hyperparameters $N'_{x1...xn}$ , there is no particular good candidate for N'
- BDeu score is very sensitive with respect to the equivalent sample size N', so need to try with many different setings
Implementation¶
let $\alpha_i = \frac{N'}{q_{i}}$¶
let $\beta_i = \frac{N'}{r_{i}q_{i}}$¶
$g_{BDeu}(G:D) = log(p(G))+\sum^n_{i=1}[\sum^{q_i}_{j=1}[log(\frac{\Gamma(\alpha_i)}{\Gamma(N_{ij}+\alpha_i)})+\sum^{r_i}_{k=1}log(\frac{\Gamma(N_{ijk}+\beta_i)}{\Gamma(\beta_i)})]]$¶
$g_{BDeu}(G:D) = log(p(G))+\sum^n_{i=1}[\sum^{q_i}_{j=1}[log(\Gamma(\alpha_i)) - log(\Gamma(N_{ij}+\alpha_i)) + \sum^{r_i}_{k=1}log(\Gamma(N_{ijk}+\beta_i)) - log(\Gamma(\beta_i))]]$¶
$g_{BDeu}(G:D) = log(p(G))+\sum^n_{i=1}[\sum^{q_i}_{j=1}[- log(\Gamma(N_{ij}+\alpha_i)) + \sum^{r_i}_{k=1}log(\Gamma(N_{ijk}+\beta_i)) )] + {q_i} \times log(\Gamma(\alpha_i)) - {r_i}{q_i} \times log(\Gamma(\beta_i)]$¶
In [9]:
class BDeuScore(BaseEstimator):
def __init__(self, data, equivalent_sample_size=10, **kwargs):
self.equivalent_sample_size = equivalent_sample_size
super(BDeuScore, self).__init__(data, **kwargs)
def BDeu_score_single(self, variable, parents, is_print=False):
var_states = self.state_names[variable]
var_cardinality = len(var_states) # number of element in a (undefined set) axiomatic set theory
state_counts = self.state_counts(variable, parents)
num_parents_states = float(state_counts.shape[1])
counts = np.asarray(state_counts)
if is_print:
printmd(f'r_i = :{var_cardinality}')
printmd(f'q_i = :{num_parents_states}')
printmd('state_counts($N_{ijk}$):'+str(counts))
display(state_counts.head(2))
print('-'*100)
log_gamma_counts = np.zeros_like(counts, dtype=np.float_)
alpha = self.equivalent_sample_size / num_parents_states
beta = self.equivalent_sample_size / counts.size
if is_print:
printmd("##### Let $alpha = {N'}/q_{i}$ = " + str(alpha))
printmd("##### Let $ beta = {N'}/(r_{i}q_{i})$ = " + str(beta))
print('-'*100)
# Compute log(gamma(counts + beta))
gammaln(counts + beta, out=log_gamma_counts)
if is_print:
printmd("$log(\Gamma(N_{ijk} + beta))$ = ")
print(f"log_gamma_counts = \n {log_gamma_counts}")
print('-'*100)
# Compute the log-gamma conditional sample size
log_gamma_conds = np.sum(counts, axis=0, dtype=np.float_)
gammaln(log_gamma_conds + alpha, out=log_gamma_conds)
if is_print:
printmd("$log(\Gamma(N_{ij} + alpha))$ = $ log(\Gamma(N_{ij} + k))$")
print(f"log_gamma_conds = \n {log_gamma_conds}")
print('-'*100)
score = (
np.sum(log_gamma_counts)
- np.sum(log_gamma_conds)
+ num_parents_states * lgamma(alpha)
- counts.size * lgamma(beta)
)
if is_print:
print(f"np.sum(log_gamma_counts) - np.sum(log_gamma_conds) + num_parents_states * lgamma(alpha) - counts.size * lgamma(beta):")
return score
def score(self, model,is_print=False):
score = 0
for node in model.nodes():
if is_print:
printmd(f"## node: {node}")
printmd(f"### parent: {list(model.predecessors(node))}")
printmd(f"score: {self.BDeu_score_single(node, model.predecessors(node),is_print=is_print)}")
print(f"="*100)
score += self.BDeu_score_single(node, model.predecessors(node))
return score
In [10]:
# Importing data
simple_df = pd.read_csv('/etlstage/PEE_joint/NUS_modules/CS5340/group_data/small.csv')
display(simple_df.head())
# making a BN
BN_1 = BayesianModel([['fare','age'], ['age','passengerclass'],['passengerclass','portembarked'], ['portembarked','survived']])
BN_1_G = nx.draw(BN_1, with_labels=True,node_size=3e3,width=3)
plt.show()
BDeuScore(simple_df, equivalent_sample_size=10).score(BN_1,is_print=True)
Out[10]:
In [11]:
# making a BN
BN_2 = BayesianModel([['fare','passengerclass'], ['passengerclass','survived'],['age','survived'], ['portembarked','survived']])
nx.draw(BN_2, with_labels=True,node_size=3e3,width=3)
plt.show()
BDeuScore(simple_df, equivalent_sample_size=10).score(BN_1,is_print=True)
Out[11]:
In [12]:
display(GME_stock_df.head(5))
# making a BN
BN_GMS = BayesianModel([['open_price','adjclose_price'],
['high_price','volume'],
['low_price','volume'],
['open_price','high_price'],
['open_price','low_price'],])
nx.draw(BN_GMS, with_labels=True,node_size=3e3,width=2)
plt.show()
BDeuScore(GME_stock_df, equivalent_sample_size=20).score(BN_GMS,is_print=False)
Out[12]:
Information-theoretic scoring functions¶
1. The score of a Bayesian network B is related to the compression that can be achieved over the data D with an optimal code induced by B.¶
2. Shannon source coding theorem. (ref: https://en.wikipedia.org/wiki/Shannon%27s_source_coding_theorem#:~:text=In%20information%20theory%2C%20Shannon's%20source,meaning%20of%20the%20Shannon%20entropy.&text=However%20it%20is%20possible%20to,with%20negligible%20probability%20of%20loss.)¶
3. Information content of D by B:¶
a. The size of an optimal code, induced by the distribution B, when encoding D¶
b. This value can be used to score the BN B. (the best BN is can be compressed close to shannon limit yet not lose info)¶
c. Gibb’s inequality¶
- Given fixed DAG struct of BN B, L(D|B)min achieved when: $\theta{ijk}=\frac{N{ijk}}{N{ij}}$
- L(T|B) is minimal when Likeihood P(D|B) is maximal
- parameters of B induces a code compresses D the most is precisely the parameters that maximizes the probability of observing D (maximum likelihood).
4. log-likelihood (LL) score function:(B->G)¶
a. This phenomenon of overfitting is usually avoided in two different ways:¶
* limiting the number of parents per network variable.
* penalization factor over LL score:
* **MDL/BIC**
* **AIC**
* **NML**
4. log-likelihood (LL) score function with penalization:¶
where $C(G)=\sum^n_{i=1}(r_i-1)q_i$ (network complexity)¶
idea:¶
ref: Journal of Machine Learning Research - SCORING BAYESIAN NETWORKS USING MUTUAL INFORMATION AND INDEPENDENCE TESTS (2154~2155)
- the best model is the one that minimizes the sum of the description length of the model and the description length of the data given the model
- For simple models require shorter description lengths, but the description length of the data given the model increases
- The minimum description length principle establishes an appropriate trade-off between complexity and precision
- for BN: the description length includes the length required to represent the network plus the length necessary to represent the data given the network
- we must store its probability values, and this requires a length which is proportional to the number of free parameters of the factorized joint probability distribution and This number, called network complexity and denoted as C(G) ## =========================================================================================
1. MDL/BIC¶
BIC: Bayesian Information Criterion / MDL: Minimal Descriptive Length¶
$f(N) = \frac{1}{2}log(N)$¶
In [13]:
class BicScore(BaseEstimator):
def __init__(self, data, **kwargs):
super(BicScore, self).__init__(data, **kwargs)
def Bic_score_single(self, variable, parents, is_print=False):
var_states = self.state_names[variable]
var_cardinality = len(var_states) # number of element in a (undefined set) axiomatic set theory
state_counts = self.state_counts(variable, parents)
num_parents_states = float(state_counts.shape[1])
counts = np.asarray(state_counts)
sample_size = len(self.data)
if is_print:
printmd(f'r_i = :{var_cardinality}')
printmd(f'q_i = :{num_parents_states}')
printmd('state_counts($N_{ijk}$):'+str(counts))
printmd('sample_size($N$):'+str(sample_size))
display(state_counts.head(2))
print('-'*100)
log_likelihoods = np.zeros_like(counts, dtype=np.float_)
# Compute the log-counts
np.log(counts, out=log_likelihoods, where=counts > 0)
if is_print:
printmd("$log(N_{ijk})(>0)$ = ")
print(f"{log_likelihoods}")
print('-'*100)
# Compute the log-conditional sample size
log_conditionals = np.sum(counts, axis=0, dtype=np.float_)
np.log(log_conditionals, out=log_conditionals, where=log_conditionals > 0)
if is_print:
printmd("$log(N_{ij})(>0)$")
print(f"{log_conditionals}")
print('-'*100)
# Compute the log-likelihoods
log_likelihoods -= log_conditionals
log_likelihoods *= counts
if is_print:
printmd("log_likelihoods = $N_{ijk} * log(N_{ijk}/N_{ij})$ = $N_{ijk}*[log(N_{ijk}) - log(N_{ij})]$ = ")
print(f"{log_likelihoods}")
print('-'*100)
score = np.sum(log_likelihoods)
score -= 0.5 * log(sample_size) * num_parents_states * (var_cardinality - 1)
if is_print:
printmd("bic score = $\sum log\_likelihoods - 0.5*log(sample\_size(N)) * num_parents\_states(q_i)*(var\_cardinality(r_i) - 1)$")
print(f"{log_likelihoods}")
print('-'*100)
return score
def score(self, model,is_print=False):
score = 0
for node in model.nodes():
if is_print:
printmd(f"## node: {node}")
printmd(f"### parent: {list(model.predecessors(node))}")
printmd(f"score: {self.Bic_score_single(node, model.predecessors(node),is_print=is_print)}")
print(f"="*100)
score += self.Bic_score_single(node, model.predecessors(node))
return score
In [14]:
# Importing data
simple_df = pd.read_csv('/etlstage/PEE_joint/NUS_modules/CS5340/group_data/small.csv')
display(simple_df.head())
# making a BN
BN_1 = BayesianModel([['fare','passengerclass'], ['age','passengerclass'],['portembarked','passengerclass'], ['passengerclass','survived']])
BN_1_G = nx.draw(BN_1, with_labels=True,node_size=3e3,width=3)
plt.show()
BicScore(simple_df).score(BN_1,is_print=True)
Out[14]:
In [15]:
# making a BN
BN_2 = BayesianModel([['fare','passengerclass'], ['passengerclass','survived'],['age','survived'], ['portembarked','survived']])
nx.draw(BN_2, with_labels=True,node_size=3e3,width=3)
plt.show()
BicScore(simple_df).score(BN_2,is_print=True)
Out[15]:
In [16]:
display(GME_stock_df.head(5))
# making a BN
BN_GMS = BayesianModel([['open_price','adjclose_price'],
['high_price','volume'],
['low_price','volume'],
['open_price','high_price'],
['open_price','low_price'],])
nx.draw(BN_GMS, with_labels=True,node_size=3e3,width=2)
plt.show()
BicScore(GME_stock_df).score(BN_GMS,is_print=False)
Out[16]:
In [17]:
class AicScore(BaseEstimator):
def __init__(self, data, **kwargs):
super(AicScore, self).__init__(data, **kwargs)
def Aic_score_single(self, variable, parents, is_print=False):
var_states = self.state_names[variable]
var_cardinality = len(var_states) # number of element in a (undefined set) axiomatic set theory
state_counts = self.state_counts(variable, parents)
num_parents_states = float(state_counts.shape[1])
counts = np.asarray(state_counts)
sample_size = len(self.data)
if is_print:
printmd(f'r_i = :{var_cardinality}')
printmd(f'q_i = :{num_parents_states}')
printmd('state_counts($N_{ijk}$):'+str(counts))
printmd('sample_size($N$):'+str(sample_size))
display(state_counts.head(2))
print('-'*100)
log_likelihoods = np.zeros_like(counts, dtype=np.float_)
# Compute the log-counts
np.log(counts, out=log_likelihoods, where=counts > 0)
if is_print:
printmd("$log(N_{ijk})(>0)$ = ")
print(f"{log_likelihoods}")
print('-'*100)
# Compute the log-conditional sample size
log_conditionals = np.sum(counts, axis=0, dtype=np.float_)
np.log(log_conditionals, out=log_conditionals, where=log_conditionals > 0)
if is_print:
printmd("$log(N_{ij})(>0)$")
print(f"{log_conditionals}")
print('-'*100)
# Compute the log-likelihoods
log_likelihoods -= log_conditionals
log_likelihoods *= counts
if is_print:
printmd("log_likelihoods = $N_{ijk} * log(N_{ijk}/N_{ij})$ = $N_{ijk}*[log(N_{ijk}) - log(N_{ij})]$ = ")
print(f"{log_likelihoods}")
print('-'*100)
score = np.sum(log_likelihoods)
score -= num_parents_states * (var_cardinality - 1)
if is_print:
printmd("bic score = $\sum log\_likelihoods - num_parents\_states(q_i)*(var\_cardinality(r_i) - 1)$")
print(f"{log_likelihoods}")
print('-'*100)
return score
def score(self, model,is_print=False):
score = 0
for node in model.nodes():
if is_print:
printmd(f"## node: {node}")
printmd(f"### parent: {list(model.predecessors(node))}")
printmd(f"score: {self.Aic_score_single(node, model.predecessors(node),is_print=is_print)}")
print(f"="*100)
score += self.Aic_score_single(node, model.predecessors(node))
return score
In [18]:
display(GME_stock_df.head(5))
# making a BN
BN_GMS = BayesianModel([['open_price','adjclose_price'],
['high_price','volume'],
['low_price','volume'],
['open_price','high_price'],
['open_price','low_price'],])
nx.draw(BN_GMS, with_labels=True,node_size=3e3,width=2)
plt.show()
AicScore(GME_stock_df).score(BN_GMS,is_print=False)
Out[18]:
2. NML scoring function (not found any implementation online, this implementation in python may be the first)¶
NML: normalized minimum likelihood¶
log-likelihood (LL) score function with penalization (General form):¶
Idea: What if we treat $C(G)f(N)$ as a distribution H?¶
- Given a set of probability distributions H the encoder thinks that there is one distribution H0 ∈ H that will assign high likelihood (low code length) to the incoming data D of fixed size N. $(L(H)_{min} = L(H0))$
- We woukd like to design a code that for all D will compress D as close as possible to the code associated to H0 ∈ H that maximizes the likelihood of D $LL_D(G)_{max} | H0$
- We call to this H0 ∈ H the best-fitting hypothesis. ## Too much explaination (ignore) ref: http://citeseerx.ist.psu.edu/viewdoc/download;jsessionid=72FFB108DFC6009541FF8E1A83E88308?doi=10.1.1.140.1665&rep=rep1&type=pdf -> implementation ## Final form of NML function: ## ${NML}(D|G) = LL_D(G) - C_N(B_G)$ where $B_G$ denote the set of all BNs with network structure G
- There is no hope for computing $C_N(B_G)$ efficiently:
- It involves an exponential sum over all possible data of size N.
- It is not decomposable over the network structure.
#### But from ref: proposed to approximate $C_N(B_G)$ by considering only the contribution to the parametric complexity of the multinomial distributions associated to each variable given a parent configuration:
In [19]:
class fNMLScore(BaseEstimator):
def __init__(self, data, **kwargs):
super(fNMLScore, self).__init__(data, **kwargs)
def Szpankowski_approximation(self, row):
if row.n == 0: return 1
return 0.5*row.n*np.pi*np.exp(np.sqrt(8/(9*row.n*np.pi))+((3*np.pi-16)/(36*row.n*np.pi)))
def Compute_regret(self, row, r):
if row.n == 0: return 1
return row[f'C_{r-1}'] + (row.n/(r-2))*row[f'C_{r-2}']
def fNML_score_single(self, variable, parents, is_print=False):
var_states = self.state_names[variable]
var_cardinality = len(var_states) # number of element in a (undefined set) axiomatic set theory
state_counts = self.state_counts(variable, parents)
num_parents_states = float(state_counts.shape[1])
counts = np.asarray(state_counts)
sample_size = len(self.data)
if is_print:
printmd(f'r_i = :{var_cardinality}')
printmd(f'q_i = :{num_parents_states}')
printmd('state_counts($N_{ijk}$):<br/>'+str(counts))
printmd('sample_size($N$):'+str(sample_size))
display(state_counts.head(2))
print('-'*100)
log_likelihoods = np.zeros_like(counts, dtype=np.float_)
# Compute the log-counts
np.log(counts, out=log_likelihoods, where=counts > 0)
if is_print:
printmd("$log(N_{ijk})(>0)$ = <br/>"+f"{log_likelihoods}")
print('-'*100)
# Compute the log-conditional sample size
log_conditionals = np.sum(counts, axis=0, dtype=np.float_)
np.log(log_conditionals, out=log_conditionals, where=log_conditionals > 0)
if is_print:
printmd("$log(N_{ij})(>0)$<br/>"+f"{log_conditionals}")
print('-'*100)
# Compute the log-likelihoods
log_likelihoods -= log_conditionals
log_likelihoods *= counts
if is_print:
printmd("log_likelihoods = $N_{ijk} * log(N_{ijk}/N_{ij})$ = $N_{ijk}*[log(N_{ijk}) - log(N_{ij})]$ = ")
print(f"{log_likelihoods}")
print('-'*100)
# Compute the 𝐶𝑁(𝐵𝐺) table
N_ij = np.sum(counts, axis=0, dtype=np.float_)
N_ij_unique, N_ij_cnt = np.unique(N_ij, return_counts=True)
N_ij_dist = dict(zip(N_ij_unique, N_ij_cnt))
C_ini = [[n, 1] for n in N_ij_unique]
C_lookup_table = pd.DataFrame(C_ini, columns = ['n', 'C_1'])
if var_cardinality > 1:
C_lookup_table['C_2'] = C_lookup_table.apply(self.Szpankowski_approximation,axis=1)
if var_cardinality > 2:
for r_i in range(3,var_cardinality+1):
Compute_regret_par = partial(self.Compute_regret, r=r_i)
C_lookup_table[f'C_{r_i}'] = C_lookup_table.apply(Compute_regret_par,axis=1)
if r_i > 5: #only obtain the last 3 regret C for memory saving
C_lookup_table = C_lookup_table.drop(columns=[f'C_{r_i-3}'])
#C_lookup_table['N_ij_cnt'] = N_ij_cnt
C_lookup_table['fNML_regret'] = np.log(C_lookup_table[f'C_{var_cardinality}'])
fNML_regret = C_lookup_table.fNML_regret.values
if is_print:
print(f"N_ij distribution:\n {N_ij_dist}")
display(C_lookup_table)
print("fNML_regret = \n"+f"{fNML_regret}")
print('-'*100)
score = np.sum(log_likelihoods)
score -= np.sum(fNML_regret)
if is_print:
printmd("fNML score = $\sum log\_likelihoods - \sum logC^{r_i}_{N_{ij}}$")
print('-'*100)
return score
def score(self, model,is_print=False):
score = 0
for node in model.nodes():
if is_print:
printmd(f"## node: {node}")
printmd(f"### parent: {list(model.predecessors(node))}")
printmd(f"score: {self.fNML_score_single(node, model.predecessors(node),is_print=is_print)}")
print(f"="*100)
score += self.fNML_score_single(node, model.predecessors(node))
return score
In [20]:
# Importing data
simple_df = pd.read_csv('/etlstage/PEE_joint/NUS_modules/CS5340/group_data/small.csv')
display(simple_df.head())
# making a BN
BN_1 = BayesianModel([['fare','passengerclass'], ['age','passengerclass'],['portembarked','passengerclass'], ['passengerclass','survived']])
BN_1_G = nx.draw(BN_1, with_labels=True,node_size=3e3,width=3)
plt.show()
fNMLScore(simple_df).score(BN_1,is_print=True)
Out[20]:
In [21]:
display(GME_stock_df.head(5))
# making a BN
BN_GMS = BayesianModel([['open_price','adjclose_price'],
['high_price','volume'],
['low_price','volume'],
['open_price','high_price'],
['open_price','low_price'],])
nx.draw(BN_GMS, with_labels=True,node_size=3e3,width=2)
plt.show()
fNMLScore(GME_stock_df).score(BN_GMS,is_print=False)
Out[21]:
In [22]:
display(GME_stock_df.shape)
display(GME_stock_df.head(5))
# making a BN
BN_GMS = BayesianModel([['open_price','adjclose_price'],
['high_price','volume'],
['low_price','volume'],
['open_price','high_price'],
['open_price','low_price'],])
nx.draw(BN_GMS, with_labels=True,node_size=3e3,width=2)
plt.show()
In [23]:
sample_size_ls = [n for n in np.arange(10, 100, 25)] + [n for n in np.arange(100, len(GME_stock_df), 100)]
print(sample_size_ls)
score_rank = {'sample_size': sample_size_ls,
'K2Score':[],
'BDeuScore':[],
'BicScore':[],
'AicScore':[],
'fNMLScore':[],
}
for sample_size in sample_size_ls:
sample_data = GME_stock_df.sample(n=sample_size, random_state=1)
score_rank['K2Score'].append(K2Score(sample_data).score(BN_GMS,is_print=False))
score_rank['BDeuScore'].append(BDeuScore(sample_data,
equivalent_sample_size=sample_size/10).score(BN_GMS,is_print=False))
score_rank['BicScore'].append(BicScore(sample_data).score(BN_GMS,is_print=False))
score_rank['AicScore'].append(AicScore(sample_data).score(BN_GMS,is_print=False))
score_rank['fNMLScore'].append(fNMLScore(sample_data).score(BN_GMS,is_print=False))
In [24]:
score_rank_df = pd.DataFrame(score_rank)
display(score_rank_df.head())
df_melt = score_rank_df.melt(id_vars='sample_size', value_vars=score_rank_df.columns.to_list()[1:])
df_melt['log_score'] = np.log(-df_melt.value)
fig = px.line(df_melt, x='sample_size' , y='log_score', color='variable')
fig.show()
In [25]:
sample_size_ls = [n for n in np.arange(10, 100, 25)] + [n for n in np.arange(100, len(GME_stock_df), 100)]
print(sample_size_ls)
score_rank = {'sample_size': sample_size_ls,
'K2Score':[],
'BDeuScore':[],
'BicScore':[],
'AicScore':[],
'fNMLScore':[],
}
c=1
for sample_size in sample_size_ls:
sample_data = GME_stock_df.sample(n=sample_size, random_state=1)
score_rank['K2Score'].append(K2Score(sample_data).score(BN_GMS,is_print=False))
if sample_size < 200:
r = sample_size/5
else:
c+=1
r = sample_size/(50)*c
score_rank['BDeuScore'].append(BDeuScore(sample_data,
equivalent_sample_size=r).score(BN_GMS,is_print=False))
score_rank['BicScore'].append(BicScore(sample_data).score(BN_GMS,is_print=False))
score_rank['AicScore'].append(AicScore(sample_data).score(BN_GMS,is_print=False))
score_rank['fNMLScore'].append(fNMLScore(sample_data).score(BN_GMS,is_print=False))
In [26]:
score_rank_df = pd.DataFrame(score_rank)
display(score_rank_df.head())
df_melt = score_rank_df.melt(id_vars='sample_size', value_vars=score_rank_df.columns.to_list()[1:])
df_melt['log_score'] = np.log(-df_melt.value)
fig = px.line(df_melt, x='sample_size' , y='log_score', color='variable')
fig.show()
In [ ]: